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ABSTRACT 

The dynamical evolution of nearly half of the known extrasolar planets in multiple-planet systems 
may be dominated by secular perturbations. The commonly high eccentricities of the planetary orbits 
calls into question the utility of the traditional Laplace-Lagrange (LL) secular theory in analyses 
of the motion. We analytically generalize this theory to fourth-order in the eccentricities, compare 
the result with the second-order theory and octupole-level theory, and apply these theories to the 
likely secularly-dominated HD 12661, HD 168443, HD 38529 and v And multi-planet systems. The 
fourth-order scheme yields a multiply-branched criterion for maintaining apsidal libration, and implies 
that the apsidal rate of a small body is a function of its initial eccentricity, dependencies which are 
absent from the traditional theory. Numerical results indicate that the primary difference the second 
and fourth-order theories reveal is an alteration in secular periodicities, and to a smaller extent 
amplitudes of the planetary eccentricity variation. Comparison with numerical integrations indicates 
that the improvement afforded by the fourth-order theory over the second-order theory sometimes 
dwarfs the improvement needed to reproduce the actual dynamical evolution. We conclude that LL 
secular theory, to any order, generally represents a poor barometer for predicting secular dynamics 
in extrasolar planetary systems, but does embody a useful tool for extracting an accurate long-term 
dynamical description of systems with small bodies and/or near-circular orbits. 

Subject headings: planets and satellites: individual (HD 12661,HD 168443,HD 190360,HD 38529,HIP 
14810) — planets and satellites: general — methods: analytical — celestial me- 
chanics 



1. INTRODUCTION 

Dynamical systems of three or more bodies typically 
include physical phenomena known as mean motion res- 
onances and secular perturbations. The former occur 
when pairs of bodies have orbital periods whose ratio 
can be approximately expressed as a ratio of two small 
integers. Otherwise, secular perturbations dominate the 
system's long-term evolution. Each phenomenon can be 
described by particular terms in a gravitational poten- 
tial which has become commonly k nown as a "disturb- 
ing function" ([Ellis fe M urray 2000), a dyna mically rich 
segm ent of the Hamiltonian of the system (Morbidclli 
2002). Laplac e-Lagrange (LL) secular th eory, reviewed 
for example in Murray fc Dermottl (|1999f ) , represents an 
elegant formulation in which the time evolution of ob- 
jects' orbital parameters can be described analytically 
by utilizing a select few terms in the disturbing function 
up to second-order in eccentricities and inclinations. 

However, the theory has come under recent scrutiny 
for its failure to reproduce the phase space struc- 
ture of dynamical syste ms of interest, such as extraso - 
lar planetary systems. iLibert fc Henrardl (|2005l 12006ft 
showcase the limitations of the theory by using an 
alternative, 12th-order exp a nsion in eccentricity, and 
iBeauge. Nesvornv. fc Donea (|2006ft successfully repro- 
duce motions of irregular satellites with eccentricities as 
high as 0.7 by using a third- order Hori averaging method. 
iBarnes fc Greenbergl j2006) find significant discrepancies 
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between first-order secular theory and N-body simula- 
tions, and lCuk fc Burns! (|2004f ) demonstrate that special 
care must be given to model correctly the se c ular m otion 
of irregular satellites. iChristou fc Murravl (|1997ft gen- 
eralize Laplace-Lagrange secular theory to second-order 
in the masses, and apply the result to the secularly- 
dominated Uranian satellite system. Several extrasolar 
systems exhibit "hierarchical" behavior, and "octupole- 
level" secular theory has b een developed in order to de- 
scribe analogous systems jK rvmolowski & Mazch 119991 : 
iFord. Kozinskv fc Rasiol [20001 : iLee fc Peald 12003ft . The 
secular phase space of multi-planet systems is complex 
(|Feiozll2002t iMichtchenko fc Malhotral I2004D . and their 
accurate characterization necessitates a comprehensive 
general analytical treatment. 

Despite thes e limitations, cl a ssical LL theo r y re- 
mains useful ([Greenbergl 119771: iHep penhcimcr 19801: 
I Chiang fc Murravl |2002t " Tji et al\ [20031: IWu fc Murra 
20031: IZhou fc Sunl 120031 IZamaska fc Tremaine! 1200 
I Adams fc Laughlinl 120061 : iNamouni fc Zh ou 2006) in or- 
der to, qualitatively at least, describe secular motion. 
The formulation of the theory admits some compact an- 
alytical results dMalhotral l2002t iMoriwaki fc Nakagawal 
120041 : IVeras fc Armitagdl2004ft and allows one to circum- 
vent singularities present in Lagrange's planetary equa- 
tions, which explicitly describe the time evolution of or- 
bital elements of a secondary. Here we extend the planar 
theory to fourth order, and demonstrate the physical con- 
sequences of this extension. In doing so, we help remove 
a major assumption (that requiring small eccentricities) 
of the standard theory, and address the limits of the the- 
ory's applicability, even in its more general state. 

At least 20 multi-planet extrasolar systems around 
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Solar- type stars are known to exist (|Butler et aZ.1 12006) . 
four of which are thrcc-planct systems (v And, HD 37124, 
HD 69830, GJ 876) and two of which are four-planet 
systems (/x Ara, 55 Cnc). The orbital periods for sev- 
eral pairs of the planets in these ps 20 systems are in a 
high-order mean motion commensurability, the effects of 
which typically dominate the long-term motion. How- 
ever, other pairs of planets are near no such commen- 
surability, and hence are described by orbital elements 
which vary primarily according to secular perturbations. 

Section 2 briefly reviews standard LL secular theory 
and introduces some notation. Section 3 generalizes the 
theory for two resonant bodies, and presents some an- 
alytical results in specific cases. The idea of "apsidal 
libration" is discussed in Section 4, and criteria for libra- 
tion are presented. Section 5 generalizes the theory to N 
resonant bodies, and Section 6 compares the second and 
fourth-order theories and octupole-level theory with nu- 
merical integrations of several extrasolar planetary sys- 
tems. We discuss the implications and conclusions of our 
model in Sections 7 and 8, respectively. 

2. STANDARD LAPLACE-LAGRANGE SECULAR THEORY 

The fully general disturbing function, Kj, where j = 
1(2) for the outer (inner) secondary, can be written as: 



hj = ej sin Wj 
ki = e.,- cos Wi 



h k = e k sin zu k 
k k = e k costUfc. 



(5) 



The disturbing function is then re-expressed in terms 
of hj and kj only, and after some algebra, one obtains 
the equations of motion: 



hi \ ( A n A 12 \ f ki 
ho / I A 21 A 22 J I k 2 



(G) 



The Laplace-Lagrange secular solution hence reduces 
to two eigenvalues and eigenvectors, which may be solved 
for exactly given the initial conditions. 

3. FOURTH-ORDER LAPLACE-LAGRANGE SECULAR 
THEORY 

In order to obtain the planar fourth-order solution, 
we retain all eccentricity terms up to fourth order in 
Eq. |l}. Therefore, we keep three arguments (with 
^(3) = 2vji — 2zu 2 ), some of which contain more than 
one relevant term. The particular functions of a = a 2 /a\ 
which the Cj represent ca n be found in Appendix B of 
iMurrav fc Dermottl (|1999l ). Using their / notation, one 
obtains: 



K, 



E 

,M=i 



i0« 



(1) 



such that i labels each secular and resonant argument 
4>, u labels each coefficient of a given argument, X( hu ) 
is a function of each secondary's inclination and eccen- 
tricity (= ei,e 2 ) alone, and cj'* u ' is a function of each 
secondary's semimajor axis (= a\,a 2 ) alone. 

In traditional Laplace-Lagrange theory, only four ar- 
guments are retained, none featuring coupling between 
eccentricity and inclination. As this work is only con- 
cerned with the former, we present only the planar theory 
in this section. The two eccentricity arguments retained 



are 



and 



vj\ — m 2 , where zu represents the 



longitude of pericenter. Manipulation yields a disturbing 
function of the form: 



K, 



Aj k e\e 2 cos (zv± — zu 2 ), 



(2) 



where the constants in time, Ajj and Aj k , are functions 
of each secondary's semimajor axis and mass (= m\ , m 2 ) , 
such that k = 1, 2 but k ^ j. This form of the disturbing 
function is then inserted into a truncated version of two 
of Lagrange's Planetary Equations (Brouwer & Clemence 
1961): 



1 



dKi 



22 dz 
Hj dj ej 

1 dK 



dej 



(3) 
(4) 



where jij = Q (m + rrij), Q represents the gravitational 
constant and m the central mass. In order to remove the 
eccentricity singularity in the denominator, the following 
auxiliary variables are defined: 



K\ =a 1 1 Qm 2 



fi + he\ + he\e\ + f 6 ef 



+/ioeie 2 cos (zux - w 2 ) 
+/i2e?e 2 cos (zoi - zu 2 ) 

+fne\el cos [2 (ru 1 - zu 2 )} 



(7) 



K 2 =a 1 1 Gm 1 



+/ioeie 2 cos (zui - w 2 ) 
+f\\e\e 2 cos (zui - zu 2 ) 

+fi7ejel cos [2 (zv\ - zu 2 )} 



We also have: 



hi = 



3 ■ h ■ 

— ej + kjZUj, 



hj hjtbj, 

e 7 " 



- a i fj - 



ej dzuj ' 



1 ~ e j dKj 
e 7 - dej 



(8) 

(9) 
(10) 

(11) 
(12) 



In the traditional Laplace-Lagrange solution, the square 
root in the numerator of Eqs. pT|) - P^)) is expanded only 
to zeroth-order. The traditional LL theory demonstrates 
that the partial derivatives should be expressed in terms 
of hj , h k , kj and k k . We can still express the more general 
disturbing function in terms of these variables only, while 
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not retaining any eccentricity terms in the denominator, 
by using the following forms: 



»i Qm 2 



2/ 2 (hi 

2j I; 1 

2N2 



+ 2/ 5 (hf + k{) (hi + k 2 2 ) 



+4/ 6 (hi + ki) + f w {h x h 2 + hk 2 ) 

+3/12 (hi + k\) (h x h 2 + hk 2 ) 
+2f 17 ^2{h 1 h 2 + k 1 k 2 f - 

(hi + ki) (h 2 2 + k\) J , 



(13) 



dK 



e 2 -^ — = a ± Qm x 



de 



2/2 (h 



2 ' k 2 ) 



+2/ s (hi + ki) (h 2 2 + hi) 



-4/ 4 (h 2 2 



+ fio (h\h 2 + kik 2 ) 



dvj\ 



dK 2 
dw 2 



+3/n (h 2 2 + k 2 2 ) (hih 2 + hk 2 ) 
+2f 17 ^2(h 1 h 2 + k 1 k 2 ) 2 - 

(hi + ki) (h 2 2 + k 2 ) 

a^ 1 Qm 2 - fi (hik 2 - h 2 ki) 
-/12 (hik 2 - h 2 k 1 ) (hi + ki) 
-4/17 (hik 2 - h 2 ki) [hxh 2 + hk 2 ) 

a{ l Qm x fw(h\k 2 - h 2 ki) 
+f 11 (h 1 k 2 -h 2 k 1 )(h 2 + k 2 ) 
+4/i7 (h\k 2 - h 2 ki) (hih 2 + kik 2 ) 



(14) 



(15) 



(16) 



Equations (f9|)- (|16jl may be manipulated in order to 
generate the 4th-order analog of the LL theory's differ- 
ential equations of motion: 



dk, + C 2 k k + kjh) 



+C^h\ + C 5 k 3 k 2 + Ctf'rfkk + C\ 3) k k h 2 



(17) 



dh 3 + C 2 h k + h 3 k 2 
+C 4 hjkl + C b h 3 h 2 + C { / ] h 2 h k + C ( 7 3) h k k 2 



+C ( s :l) h j k j k k + Cghkkjkk + C^hj 



(18) 



where all the Ci, i — 1...9, are constant functions of 
a, and the scaling factors C„ are constant functions of 



both secondary masses and semimajor axes. These con- 
stant values are provided in Appendix A. This form of 
the equations illustrates well the dependence of all four 
variables on one another, the symmetry of those vari- 
ables, and the asymmetry in some of the constants. By 
retaining only the first two terms in the brackets of these 
equations, one recovers traditional Laplace-Lagrange sec- 
ular theory. 

Equations (fl~7|) - (|18|) are a set of four first-order cou- 
pled highly nonlinear differential equations, and in gen- 
eral need to be solved numerically. However, in the spe- 
cial case where one secondary on a circular orbit is very 
massive compared to the other secondary, we can obtain 
an analytic result. Assuming e k = for all time, we need 
only to solve: 



nti) l. . 

0„ fv, 



k 3 = -C^h 3 



Ci + C { / ] (h 2 + k 2 ) 
Ci + C { 3 3) (h 2 + k 2 ) 



(19) 
(20) 



This set of equations is also nonlinear; however, the 
presence of the squared quantities inside the parenthesis 
suggests a sinusoidal solution. Assuming body j's initial 
eccentricity equals eo, the solution is: 



hj(t): 
kj(t): 



eo sin 



eo cos 



C« (a 



00 

0^3 



elC, 



(3) 
0^3 



elC 



(21) 
(22) 



and therefore body j's longitude of pericenter varies 
linearly with a constant of proportionality equal to 

(Ci + e^C^p^ . In the traditional Laplace-Lagrange 

solution, this rate is independent on the initial eccentric- 
ity. We can analyze this rate further by expressing it in 
terms of orbital elements (details of the derivation are 
found in Appendix B): 



Ci + e$C 3 



2/^(1) . 



H a 

32 

r 1 ^2^(2) _ 3 2 

45 4 

H a 

32 



-a 2 (1 + 2eg) 



45 



1 - 



(l + ^)+0(a 6 ) 



(23) 



(24) 



The above equations demonstrate that the secular per- 
turbations between a massive body on a circular or- 
bit and a much less massive eccentric exterior body 
will cause the latter 's pericenter rate to increase with a 
greater initial eccentricity. Alternatively, a less massive 
eccentric body interior to a massive object will precess 
at a slower rate given a greater initial eccentricity. Such 
dependencies are not present in the traditional Laplace- 
Lagrange secular theory, which predicts (under the as- 
sumption of a massive central body), the rates to be 
equal. 

Furthermore, by making use of the critical planetary 
separation needed to achieve Hill Stability (|Gladmanl 
1993?; IVeras fc Armitagell2004D . we can derive an upper 
bound for the rate of change at pericenter or precession of 
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a perturbed object exterior to a massive planet (orbiting 
a much more massive star), in units of radians/year: 



(^l)r 



3miOo 



4 a 2 + 2-36TO 1 3 



(25) 



where mi is measured in solar masses and a 2 is measured 
in AUs. For mi = mj up and a 2 = 1 AU, (wi) max gives 
a maximum rate corresponding to a ~ 16, 000 yr period. 

4. APSIDAL LIBRATION AMPLITUDES 

In traditional LL theory, one nonzero argument in 
the disturbing function exists whose time derivative is 
equal to (mi — T&2), a quantity known as the "apsidal 
rate". When this rate librates, rather than circulates, 
it is sometimes claimed that the system is in "apsidal 
resonance", which has been the subject of much study 
(iBreiterl 119991: iMalhotral l2002t IChiang fe Murravl 120021: 
Beauee. Ferraz-Mello fc Michtc henkoll2003t IZhou fc Sunl 
2003h . As observed bv lBarnes fc Greenbergj (|2006f ). the 
definition of "resonance" has been multiply -defined in the 
literature. Figure 9.3 of lMorbidelill (|2002f) demonstrates 
that in fact libration and resonant regions of phase space 
may be distinct and/or independent. We don't expound 
on this matter here, but instead observe that libration 
and resonance are typically closely linked, and that an 
account of the former often aids in describing the latter. 

Librational motion requires that the sign of (jr- 2 ^ 1 (= 

mi — m 2 ) change. The amplitude of libration, often used 

to indicate how "deep" the system is in resonance, may 

be found by evaluating the value of <fy- 2 > satisfying (p 2 > = 

(2) 

0. The resulting amplitude, <fi max, IS: 



COS 
Cl 

c 2 



0l 2 l) 



cP 



e\e 2 



(26) 



Equation (I26[) may be compa r ed di rectly with Eqs. 
(16)-(21) of iBarnes fc Greenbergj (f2006h . The derivation 
of Eq. (j2"B")) suggests that no apsidal libration occurs if 
either body's orbit is circularized, and even if both ob- 
jects harbor eccentric orbits, the absolute value of the 
right-hand-side of the equation must be less than unity. 
These represent the conditions for apsidal libration to 
occur according to standard LL theory. 

For equal mass bodies, 



cos 



Ci 
C 2 



^e 2 



y/a) eie 2 



(27) 



6 (2) 
■max 



whereas for equally eccentric bodies, 
cos -1 (Ci/C 2 ). One must take care to not approx- 
imate the ratio Ci/C 2 by low-order powers of a, as 
convergence can be slow, requiring several tens of terms 
of a. A numerical sampling of phase space indicates 
that this ratio is less than unity when a > 0.409. This 
result is independent of the masses or eccentricities, 
assuming that the eccentricities of both bodies are 
equal. Physically, for separations greater than those 
given by this criterion, <f)^ must circulate, as in this 



case, the bodies are too far from each other to be "in 
resonance" . On the other extreme, if the bodies are too 
close to each other, the system will become unstable 
due to the saturation of mean motion resonances. 
Wisdom (1980) provides an approximate criterion for 
this critical distance, which is proportional to the 2/7ths 
power of the primary-secondary mass ratio. If a lies 
between these two extremes, and is sufficiently far from 
any sufficiently strong mean motion resonance, then 
standard LL theory predicts apsidal libration will occur. 

In the fourth-order theory, many more terms are in- 
cluded in the disturbing function. However, the relation 
0( 3 ) = 2</>( 2 ) allows us to derive the apsidal librational 
amplitude: 



COS ('/'max) 

(c 2 + 



c 



(2) e 2 



ci x) (c - c, 



(1)„2 



2C 9 eie 2 ( CP - cf 5 



-{28) 



where 



D = 



(c 2 + CPe 2 )-CP (c 2 + CPe 2 
4C 9 (C« - C«) \e 2 2 CP (Ci + C^el + C 4 e 2 ) 



AC? (d 



C 4 e 2 ) 



(29) 



The different solution branches manifest themselves in 
a myriad of ways, and are multiply-dependent on the ec- 
centricities and semimajor axes. Recall that in the limit 
of one planet of negligible mass on an eccentric orbit, 
its eccentricity remains constant, while its longitude of 
pericenter precesses at a constant rate, thereby prevent- 
ing apsidal libration. Therefore, Eq. (|28|) only applies in 
the case of two nonzero eccentricities. Yet, the equation 
is well-behaved as a circular orbit is approached. 

5. N-BODY FORMULATION 

Extrasolar planetary systems with at least three plan- 
ets are being discovered in increasing numbers, and the 
Solar System provides an example of a system in which 
several terrestrial and giant planets may coexist in a con- 
figuration stable over timescales of at least several Myr. 
A secular theory which could accurately describe such 
systems would be useful. We here present the fourth- 
order LL theory for N bodies. The disturbing function 
for secondary j, j = 1...N, is: 



N f 



e j e k 



[S(j-P)ft ) + S(k- P )ft ) 
* + fw^ e j e k cos (mi - m p ) 
~SU~p)ft[ p) +S(k-p)ft l 2 p) 



cos (mi — m p ) + /17 e|e| cos [2 (mi — m p )] \ (30) 
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where I = min (j, k), p = max (j, k), and ff = 
fi(a p /ai). After following the same procedure used to 
derive Eq. (fT5| . one obtains, 



AT 

E 



C (i P ) 



C-^ kj ~\~ kf» H - C3 ^j^j 



-Cg p ^k 2 k k 



Cg p ^kjkl 



■ d i lv) k k h 2 j + c<i p) k j h j h k 



+c i j p) k k h ] h k + c'i p) k] 



(31) 



E c ^ 

k =i, k =£i 



C[ lp) h 



C^hk + C^h.k] 



+C^ hj kl 
+Ct ] h 2 jh k 



d lp) hih 



3 n k 



C<r lp) h k k? 



C'q hj kj kjg 



(32) 



with new constants labeled Cf p ^ whose explicit form is 
provided in Appendix A. In order to acquire a qualita- 
tive perspective on the secular motion of an eccentric 
terrestrial planet in the midst of several giant planets 
on circular orbits, assume all bodies except body j are 
massive with zero eccentricities. Then, 



hj(t) = e sin £ 



-\ r (ip)( r ( 



%(i)=e co 8 [£jl Wj .Ci W (<7< 



Op) _ 
Op) 



0p)\ 



(33) 



The criteria for libration can be evaluated using the 
N-body formulation to determine, for example, if the ap- 
sidal libration of terrestrial planets interior to a giant 
planet is greater than if the terrestrial planets were ex- 
terior to the giant planet. Consider a pair of terrestrial 
masses in the same system with a giant planet, and de- 
note M. as the terrestrial planet/giant planet mass ra- 
tio. The giant planet may reside in between (= Case 
1), exterior to (= Case 2), or interior to (= Case 3) the 
terrestrial planets. In all cases, the planets are labeled 
such that a\ > a 2 > a 3 , and the arguments <j> such that 
(f>i p = {vji — n7 P ) max - According to the second-order the- 
ory, and assuming the giant planet is much more massive 
than each terrestrial planet, the resulting amplitudes of 
the librating angles may be expressed as: 



Case 1 : 

COS 013 = 

Case 2 : 

COS 023 = 

Case 3 : 

COS 012 = 



2 

M 

2 

M 

2 

M 



( 23 ) „2 



(12) , 2 

J 2 V«i3 a i2ei 



/i ( o 3) ai2 (v 7 "!? - 1) eie 3 



r(!3) 2 



(12) , 2 



f: 



( 23 ) c 2 



l) e 2 e 3 



(13) ; 2 

J 2 \fotv2Oix2ei 



/i ( o 2) ai2 (V"12 - l) eie 2 



(34) 



(35) 



(36) 



where we have defined ai p = a p /ai. We can use Eqs. 
(j3"4")) - (131)1) to sample regions of semimajor axis and ec- 
centricity phase space in which apsidal librational may 
occur. For either an interior or exterior pair of terrestrial 
planets, libration may only occur if they are close to each 
other but far from the giant planet - in effect, when the 
terrestrial planets represent an isolated system. In no 
instance does a giant planet embedded in between the 
terrestrial planets allow for apsidal libration. 

Although such formulas may apply in Solar System 
analogs, with giant planets on near-circular orbits, 
known extrasolar systems commonly harbor giant 
planets on moderately or highly-eccentric orbits. The 
origin of these eccentricities remains unclear, but 
may be reproduced based on gra vitational scattering 
among mu l tiple m assive planets (iRasio fc Fordl 1 1996 
Linfcldal fl997t lL evison. Lissaucr & Duncan! 11995 



Marzari fc Weidenschillind l2002t iFord 



20031: iFord. Lvstad fe Rasiol 120051 : iVeras fc Armitage 
2006T) or di sk-planet interactions (lArtvmowicz et al.l 

199ll: iPapaloizou. Nelson fc Massed 120011 : 

Qgilvie fc Lubowl 120031 : iGoldreich fcljarl l2003h . If 
we assume that a pair of terrestrial planets coexists with 
an eccentric giant planet, then the angles (071—073), 
(n72 — 1773) and (vo\ — 072) are all involved in the criteria 
for libration. The maximum librational angles may be 
manipulated into the following form: 



C0S ^ = q«p) 2 + sop) 2 ' 

(37) 

where formulas for Q<- lp \ and S ilp \ which are all 

functions of ais, 023, ai2) ei, e^, and one additional 
angle, are provided in Appendix C. This additional an- 
gle is equal to (w\ — V02) for p — 3 and {vj\ — w^) for 
p = 2. The free parameter introduced by this additional 
angle is multiplied by the eccentricity of the giant planet, 
which is assumed to be zero in Eqs. (j3"4"|) - (f3"6")l . The free 
parameter has a marked effect on the subsequent phase 
space analysis, and such an analysis may encompass an 
entire study. 

6. APPLICATION TO KNOWN EXOSYSTEMS 

6.1. Overview 

Extrasolar planet data is updated constantly, and the 
resulting orbital solutions are bound to change with the 
passage of time. Nevertheless, in many cases the er- 
rors on currently known orbital parameters allow one 
to achieve a representative qualitative description of the 
motion, and often a highly accurate quantitative solu- 
tion as well. In order to help determine which extrasolar 
planets are dominated by secular perturbations instead 
of mean motion resonances (MMRs), we compute the 
"effective distance" away from a strong MMR for each 
pair of planets. The "nominal" (approximate) locations 
of MMRs are given by [(p + n)/p]~( 2 / 3 ), where n rep- 
resents the order of the resonance and p is an integer. 
Therefore, a scale-free method of representing this "ef- 
fective distan ce" from resonance is to use this quantity 
as the metric (jChampenois fc Vienndfl999l) . Figures [T]|5] 
provide a summary of our current (as of Jan. 1, 2007) un- 
derstanding of all known multiple-planet extrasolar sys- 
tems. A schematic of the effective distance from reso- 
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47 UMo HD 74156 



Fig. 1. — Effective distance from mean motion resonance for 
all known pairs of planets in two-planet multiplanet systems, mea- 
sured in semimajor axis ratio. Locations of 1st, 2nd, 3rd, 4th, 5th 
and 6th order mean motion resonances are given by dots of de- 
creasing size. The lower horizontal line is on a broader scale than 
the top line. 



HD 69830 b,c 55 CnC b,e HD 69830 c,d 




Ups And c,d 



mu Aro c,e 

Ups And b,c 
HD 69830 b,d 
GJ 876 <:,d 55 CnC c,d 



mu Ard d,e 
55 CiC c,e 

GJ 876 b,d 
HD 37124 b,c mu Aro b,d 



Fig. 2. — Effective distance from mean motion resonance for all 
known pairs of planets in multiplanet systems with more than two 
planets, measured in semimajor axis ratio. Locations of 1st, 2nd, 
3rd, 4th, 5th and 6th order mean motion resonances are given by 
dots of decreasing size. The lower horizontal line is on a broader 
scale than the top line. 
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TABLE 


1 









Summary of masses, semimajor axes, eccentricities 
and longitude of pericenters for the multi-planet 
extrasolar systems studied here. these parameters 

were taken from jean schneider's extrasolar 
Planets Encyclopedia with sini assumed to be unity. 

nance is provided for all two-planet systems (Fig. [1} and 
all known systems with more than two planets (Fig. [2]). 
Table 1 lists the relevant orbital parameters for the sys- 
tems simulated in this study. We deem these systems to 
be sufficiently far from any MMR, and hence "secular" . 
All orbital parameters are taken from Jean Schneider's 
Extrasolar Planets Encyclopedia 3 , and for all data we 
take sin i to be unity. All N-body numerical integrations 
are carried out with the Bulirsch-Stoer HNbody integra- 
tor (K.P. Rauch & D. P. Hamilton 2007, in preparation). 

We study only those planets which satisfy a particu- 
lar constraint regarding Kaula's (1961; 1962) expansion 
of the distur bing funct i on, w hich is based on Laplace 
coefficients. ISundmanl (|1916f ) determined that in the 
planar form of this expansion, Laplace coefficients are 
sure to converge only if the following criterion is satisfied 
(|Ferraz-Mellolll994[Sidlichovskv fe Nesvornvlll994f ): 



a 2 -ff(e 2 ) < a\h{ei), 



where 



H(q) = \J 1 + q 2 cosh w + q + sinhiu, 
h(q) = \J 1 + q 2 cosh w — q — sinhtt;, 



(38) 

(39) 
(40) 

such that w is implicitly defined as w — qcoshw. We ap- 
plied this test to each pair of planets in the same system; 
failing this test implies that the use of LL theory, or any 
other theory relying on Laplace coefficients, is suspect. 
Therefore, we only consider pairs of planets which pass 
this test. 

6.2. Systems with two planets 
6.2.1. HD 12661 

iLee fc Peald (120031) use o ctupo le-level secular theory 
and iRodriguez fc Gallardol (|2005t ) use a sixth-order ex- 
pansion of Kaula's disturbing function in order to study 
the planets orbiting HD 12661, and conclude that the 
system is strongly dominated by secular motions. How- 
ever, older work indicated that the system migh t be 
trapped in a 6:1 (fGozdziews ki fc Macieie wski 200 3]) res- 
onance or is on the border of a 11:2 (Gozdzie wskill200"3T) 
MMR. 

3 at http://vo.obspm.fr/exoplanetes/encyclo/encycl.html 



Generalized Secular Theory 



7 



<£ 180 

S 120 

e 60 

c 



-o -60 

•a -120 



0.401 
0.351 
0.302 
« 0.252 
0.203 

0.154 
0.104 
0.388 
0.338 







.•/•A 

' '■ n 

■ *\ 

• A : 


t ■ A > '■■ A 1 


// v\ \ 1 

1 - x \ ' ' J 


A r 

A /' 

- A/< 

- 'A . '' 


■ A ■ / 

■ 1 \ J 

■ I.' 
' \' ' 


■ A / ' '■ A \ ' 
\ h i ; \ 1 


A :' 1 \ '• / 


u 


v'w 


■■: * / • • v ' 






0.000 0.025 0.050 0.075 

Time (Myr) 



0.100 



Fig. 3. — Eccentricity and apsidal angle variation according to 
the traditional LL theory (solid lines in the lower two panels and 
leading dark loops in the uppermost panel), the fourth-order LL 
theory (dotted lines in the lower two panels and light loops in the 
uppermost panel), and octupole-level secular perturbation theory 
(dashed lines in the lower two panels and trailing dark loops in the 
uppermost panel) for the HD 12661 system. 

We seek to determine how well the traditional and 
generalized versions of LL theory reproduce motion in 
HD 12661, and further compare these methods with 
octupole-level secular theory. Figs. [3JH] explore the dy- 
namical state of the system, and demonstrate that none 
of the theories employed here model the true evolution 
well. Considering first the variation of the eccentricity 
of ei, the fourth-order LL theory matches best the sec- 
ular periodicity, eccentricity amplitude, and eccentricity 
values of the profile exhibited by the true evolution. Oc- 
tupolc theory matches fourth-order LL theory in their 
predictions for the lower bounds of the eccentricity of 
both planets, and both theories represent a significant 
improvement over traditional LL theory with regard to 
eccentricity amplitude. All theories and the actual dy- 
namical evolution exhibit apsidal libration with a ~ 60° 
amplitude. Overall, the profiles and periodicities of the 
apsidal libration for both the fourth-order LL theory 
and octupole theory improve upon the traditional the- 
ory, with the frequency from the fourth-order LL theory 
mirroring the actual evolution best but still not well. 

6.2.2. HD 168443 
Although the HD 168443 system (jMarcv et all 120011) 



is thought to harbor a brown dwarf (jUdrv et all 120021 : 

iReffert fc Quirrenbachll2006f) . lErdi et all (|2004f ) investi- 
gate the dynamical structure of a possible habitable zone 
in the system. Even though HD 168443 b has a large ec- 
centricity (e2 « 0.53), the wide separation of both plan- 
ets allows the Sundman inequality to be satisfied. This 
same eccentricity, however, produces a significant differ- 
ence in the dynamical predictions from the second- and 
fourth-order LL theories (Fig. Both the fourth-order 
LL theory and the octupole theory reproduce the eccen- 
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Fig. 4. — Eccentricity and apsidal angle variation from an N- 
body simulation for the HD 12661 system. 

tricity values and amplitude of both planets well, even 
including the short-period variations in the true evolu- 
tion of e\ (Fig. [6]). Despite this good agreement, over 
the course of just 10 5 yr, the frequency of the eccentricity 
variation drifts by about half of a period. The circulating 
apsidal angle observed in the true evolution has an ad- 
ditional small-period noisy variation which encompasses 
both the profiles from the fourth-order LL theory and 
the octupole-level theory over 10 5 yr. 

6.2.3. HD 38529 

The HD 38529 system has been studied in the same 
context as the HD 168443 system, as both contain a 
planet which could be a brown dwarf. In HD 38429, how- 
ever, the eccentricity of the more massive, outer, planet 
exceeds that of the inner planet. The result of this dif- 
ferent configuration is that octupole theory matches the 
true evolution better than fourth-order LL theory, as can 
be seen in Figs. [7]|S] The agreement achieved in the 
eccentricity values and amplitudes of the inner planet 
is particularly good (within 1%) between the octupole 
theory and the true evolution. The eccentricity ampli- 
tude achieved by the fourth-order LL theory is actually 
worse than that from the traditional theory, although 
the secular frequency is better modeled by the fourth- 
order theory. The long (> 5 x 10 4 yr) circulation period 
of the apsidal angle differs from that predicted by the 
fourth-order LL theory and octupole theory by only a 
few percent. 

Despite some positive agreement, no theory varies the 
outer planet's eccentricity by more than 0.0003, and this 
eccentricity is centered nearly around = 0.3600. Both 
these characteristics starkly display the inadequacy of 
any of these theories to quantitatively describe the ad- 
mittedly minor but non-negligible perturbations by the 
less massive inner planet on the more massive moderately 
eccentric outer one. 
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Fig. 5. — Eccentricity and apsidal angle variation according 
to the traditional LL theory (solid lines in the lower two pan- 
els and light squares in the uppermost panel), the fourth-order 
LL theory (dotted lines in the lower two panels and dots in the 
uppermost panel), and octupole-level secular perturbation theory 
(dashed lines in the lower two panels and thick triangles in the 
uppermost panel) for the HD 168443 system. 
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Fig. 7. — Eccentricity and apsidal angle variation according 
to the traditional LL theory (solid lines in the lower two pan- 
els and light squares in the uppermost panel), the fourth-order 
LL theory (dotted lines in the lower two panels and dots in the 
uppermost panel), and octupolc- level secular perturbation theory 
(dashed lines in lower two panels and thick triangles in the upper- 
most panel) for the HD 38529 system. 
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Fig. 6. — Eccentricity and apsidal angle variation from an N- 
body simulation for the HD 168443 system. 



6.3. Systems with more than two planets 

Although one pair of planets in a multi-planet system 
may reside in a mean motion resonance, other pairs of 
planets in the system may evolve primarily by secular 
means. In three-planet systems with a MMR between 
two of the planets, the other planet's evolution will be af- 
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Fig. 8. — Eccentricity and apsidal angle variation from an N- 
body simulation for the HD 38529 system. 

fected, albeit indirectly, by non-secular motions. Strong 
(1st, 2nd, or 3rd order) MMRs are likely to significantly 
influence the \i Ara, GJ 876, and 55 Cnc systems, and 
planets in the 55 Cnc, GJ 876, v And, and HD 69830 sys- 
tems are all close enough (a < 0.07 AU) to their parent 
stars to have been tidally influenced, possibly resulting 
in a reduction of their eccentricities. Further, a pair of 
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Fig. 9. — Eccentricity and apsidal angle variation according 
to the traditional LL theory (solid lines in the lower two pan- 
els and light squares in the uppermost panel), the fourth-order 
LL theory (dotted lines in the lower two panels and dots in the 
uppermost panel), and octupole-level secular perturbation theory 
(dashed lines in the lower two panels and thick triangles in the 
uppermost panel) for the planet pair v And c and v And d. 
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Fig. 10. — Eccentricity and apsidal angle variation from an 
N-body simulation for the planet pair v And c and v And d. 

planets in the HD 37124 system fail to satisfy the Sund- 
man Criterion for convergence of Laplace coefficients in 
terms of the disturbing function. All these observations 
cast doubt on the effective use of LL secular theory, even 
a generalized N-body version, on the majority of the cur- 
rently known three and four-planet extr asolar planetar y 
systems. One exception is provided by Ui et all (|2006|) . 




0.25 0.50 0.75 

Eccentricity of Outer Body 

Fig. 11. — Regions of the v And c - v And d eccentricity phase 
space which allow for libration from the traditional LL theory only 
(light gray), traditional LL t heo ry and the generalized LL theory 
for the lower sign from Eq. (1 28 I t only (medium gray), traditional 
LL t heo ry and the generalized LL theory for the upper sign from 
Eq. d 28 1 ) only (dark gray), and the traditional LL theory and both 
branches of the generalized LL theory (black). 

who demonstrate that the orbital evolution of the three 
planets in the HD 69830 are well-described with LL secu- 
lar theory due to their small eccentricities and relatively 
large orbital separations. 

Despite these reservations, one may apply LL theory 
to the planet pair v And c and v And d. The results, pre- 
sented in Figs. IMTTl perhaps showcase best how even the 
fourth-order theory can fail. Both the traditional LL the- 
ory and the fourth-order theory predict circulation of the 
apsidal angle, whereas only the octupole theory instead 
correctly predicts apsidal libration. This key distinction 
casts doubt on the ability of fourth-order LL theory to 
make reliable qualitative or quantitative predictions for 
multi-planet extrasolar systems in general. 

Still, such systems may demonstrate the variety of 
phase space portraits allowed by the theory. Figure 
1111 displays a panorama of regions in eccentricity phase 
space in which apsidal libration is allowed according to 
traditional LL theory and fourth-order LL theory. The 
figure hints at the branch point introduced by the dis- 
criminant in Eq. (|28|) , and demonstrates that traditional 
LL theory overestimates the region of phase space which 
allows for apsidal libration to occur. This dense and var- 
ied phase portrait demonstrates how strongly the fourth- 
order modification to traditional LL theory may play a 
role in an extrasolar system. 

7. DISCUSSION 

Our demonstration of LL theory's limitations 
for extrasolar planetary motions of known sys- 
tems should not overshadow its utility in modeling 
other dynamical systems, including Solar Sys- 
tem analogs. Although MMRs help establish the 
framework of planetary systems (|Gomes et ~aH [20051: 



10 



Veras & Armitage 



iMorbidelli eTal\ 120051: iTsiganis et o/]|2005[ ). their final 
states m ight be dominated by purely secular pertur- 
bations (iBrouwer fc van" Woerkoml fl950L lLaskarl [1988; 
llto fc Tanikawall2002l ). Quantitatively accurate secular 
solutions are difficult to obtain, even on a ~ 1 — 100 
Myr timescale, as both stable and chaotic long-term 
(> 1 Gyr) solutions exist for the Solar S ystem withi n 
the uncertainties of current measurements (Hayes 2006). 
Therefore, analytical results such as those from Eqs. 
f2"I j) -(|2"3 |) . (USD and dMD-flU may give re P resentative 
qualitative descriptions of the dynamical state of the 
system on timescales shorter than those which may give 
rise to destructive stochasticity. 

A purely secular theory, by definition, does not incor- 
porate the effects from the resonant terms which may 
cause the theory to fail. The true disturbing function 
includes an infinite number of terms, and depending on 
the orbital parameters of both resonant bodies, each term 
containing mean longitudes has a particular "strength" . 
In essence, a secular theory eliminates all of these terms, 
and justifies this neglect by appealing to how the com- 
bined strength of these terms should not match the 
strength of those terms without mean longitudes. In real- 
ity, however, resonant terms play a non-negligible role in 
the evolution of seemingly secular planetary and satellite 
s ystems. 

iMalhotra et ali ()1989l ) presented a way to incorpo- 
rate "near-resonant" terms into Laplace-Lagrange secu- 
lar theory, significantly improving their results. A similar 
approach could be invoked for fourth-order LL theory in 
order to improve its current general lack of agreement 
with the actual evolution of multi-planet exosystems. If 
the strong "near-resonant" term in question is a low- 
order resonance, then only one or two resonant terms 
might significantly impact the motion. With assump- 
tions about the period of these terms, their tendency to 
circulate instead of librate at near-commensurabilities, 
the constancy of their coefficients, and their decoupled 
states, one can model each term as a pendulum, the 
paradigm for single-resonance models. The pendulum 
model allows one to define an energy and compute a cir- 
culation period for each resonant term, and ultimately to 
include th e relevant avera g ed eff ects into the disturbing 
function. IMalhotra et al\ (|1989h perform this task for 
first-order near-resonances. As per their discussion, ex- 
panding their theory in order to encompass a treatment 
such as ours would involve deriving the corrections due to 
second-order near-resonances and possibly dealing with 
the corrections introduced from coupling between reso- 
nant terms. 

Because secular eigenfrequencies are proportional to 
the resonant mass/central mass ratio and to the mean 
motion of the resonant masses, the utility of LL theory 
largely hinges on the values of the masses and semimajor 
axes of the r esonant objects . Unl i ke for secular satel- 
lite systems (jMalhotra et all 119891: IChristou fc Murray! 
[1993), the mass ratio in extrasolar systems is orders of 
magnitude larger. The result is that the quality of the 
results obtained from LL theory is sensitive to this ratio, 
a result also suggested by our numerical investigations. 
The planets in HD 108873 are close to a 4:1 MMR. Us- 
ing the parameters for the planets from Table 1, we find 
that the system's properties are sensitively dependent on 
the exoplanet masses and semimajor axis ratio, but are 



more robust to changes in initial eccentricity. A possible 
(but speculative) reason for this characterization is that 
eccentricities appear in the expressions for coefficients of 
resonant terms, but are not included in the expressions 
for the secular eignefrequencies. 

The neglect of inclination in our fourth-order theory 
- although already justified by the 1) unknown inclina- 
tions of the vast majority of extrasolar planets, 2) near- 
coplanarity of several dynamical subsystems in the Solar 
System, and 3) complexities arising from the coupling of 
eccentricity and inclination in fourth-order LL theory - 
may not be of significant consequence given the marginal 
improvement afforded by expanding the theory to include 
eccentricity terms up to fourth-order. However, such 
coupled terms might introduce important nonlinearities 
in the solution of LL theory for fully three-dimensional 
problems, and deserve consideration in future studies. 

A subsequent planar expansion to arbitrary order is 
possible and may proceed along the same lines as the 
method outlined here. Expansion to Ath order will 
involve terms in the disturbing function of the form 
cos [N (wi — TUp)] , which can be expressed as polynomi- 
als in cos raj; and cosm v . The differential equations repre- 
senting the motion will contain hj, hk, kj and fcfe terms of 
all odd orders up to (N — l)th order. However, such an 
expansion may only be useful when seeking a highly ac- 
curate solution in systems where effects from barycentric 
drifts, short-period terms and MMRs may all be safely 
neglected. 

One may also generalize LL theory to fourth-order in 
the inclinations but not the eccentricities when modeling 
inclined circular orbits. Under the small angle approx- 
imation, sin/j/2 ss Ij/2, where Ij is the inclination of 
body j, the equations of motions and criteria for nodal 
libration will have exactly the same form as those of the 
planar case, except that the constants will not depend on 
the relative locations of the bodies being perturbed (i.e. 
no delta functions would appear in the analog of Eq. I3T))) 
as there is no preferred reference plane. A consequence 
is that the nodal rate of an inclined massless particle 
interior to a giant planet on a circular coplanar orbit 
would remain unchanged if the bodies switched semima- 
jor axes. To fourth order, this rate would be dependent 
on the massless particle's inclination, and be equal to: 
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(41) 



Comparison with Eq. ([24]) demonstrates a difference 
of sign, but a remarkable similarity of coefficients. In 
fact, reduction to the traditional LL theory and to second 
order in a demonstrates that the nodal rate exceeds the 
apsidal rate by exactly a factor of a. The differences in 
the presence of powers of a in Eqs. (|24|) and ([4Tjl are due 
to the different functions of the Laplace coefficients which 
arise naturally from the expansion of Kaula's disturbing 
function about zero eccentricities and inclinations. 

8. CONCLUSION 
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APPENDIX 
APPENDIX A 



The following constants in time, expressed as functions of the planetary scmimajor axes and masses only, are used 
in the two-body fourth-order LL theory from Eq. (fT%|) : 



C'p = Ga 1 2 m 2 n 1 
C x = 2/ 2 



c 4 = 

= 3/ 12 - i/ : 



= 4/ 6 - h 
2/ 5 - 2/l7 

12 — j/lO 



CP = 2/ 12 
Cg = 4/ 17 



C( a ] = Q a l lfl 2 2 m 1 H 2 
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^7 — 111 



1 _ 7^/10 
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The following constants are a modified, more general, form of those in Eqs. (|A1[) - 
order LL theory from Eq. 



used in the N-body fourth- 
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APPENDIX B 



One may express constant s found in the coefficients of terms in the disturbing function as functions of derivatives 
of Laplace coefficients. From lMurrav fc Dermottl (|1999l ). secular dynamics yield the following forms: 



h = \aDb^ (a) + la 2 V%W (a) 
42 52 

1 



h = --ab { - 1] (a 
4 2 



-abf (a) 

4 2 



32 2 LZa 2 



(Bl) 
(B2) 
(B3) 
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h = (a) + ^-a 2 V 2 b^ (a) + -a^b^ (a) + ±-a*V% { ? (a) (B4) 

2 lb 2 4 2 62 2 

f G = LcCDtf) {a ) + ^-a 2 V 2 b^ ( a ) + |-a 3 P 3 ^ (a) + ±=oWb<? (a) (B5) 
lb 2 62 2 62 2 128 2 

/ 8 = la 2 fe(- 2 > (a) + \a 2 bf (a) + ^-a 2 bf (a) (B6) 

16 2 4 2 16 2 

ho = UP (a) - (a) - \a 2 V 2 bP (a) (B7) 

Z 2 Z 2 4 2 

/n =-\a 2 V 2 b^ (a) ^-a 3 V 3 b^ (a) - (a) (B8) 

o 2 lo 2 oz 2 



12 : 



IftW ( a ) _ I a 2? 5 M (a) - ^a 2 V 2 b ( P (a) (B9) 
82 82 16 2 

- A a 3p3 6 (D (a) _ _L a 4p4 & (l) (Q) (B10) 
ID 2 2 

/ 17 = A&< 2 > (a) - A^ftPO (a) + ^a 2 ^ (a) (Bll) 
16 2 lb 2 62 2 

+ ia 3 P 3 6 (2) (a) + — a 4 P 4 6 (2) (a) (B12) 
82 64 2 

where 61 is a Laplace coefficient and I? is a differential operator. Using the hypergeometric expansion of the Laplace 
coefficients ([Murray & Dermottl[l999l) one may show that: 

„ 3 2 45 4 525 , 11025 8 218295 1(1 2081079 12 

f 2 = -a 2 H a 4 H a 6 + a 8 + a H a 12 B13 

8 64 512 8192 131072 1048576 1 ' 

3 2 45 4 525 a 11025 „ 218295 in 2081079 12 

f 3 = — a 2 a 4 a 6 a 8 a 10 a 12 B14 

J3 2 16 128 2048 32768 262144 1 ' 

, 135 4 2625 fi 231525 8 1964655 10 114459345 12 

f 4 = a 4 + a 6 H a 8 H a H a 12 B15 
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9 2 225 4 11025 fi 99225 s 12006225 in 81162081 12 
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, 15 2 945 4 4725 fi 606375 8 8513505 10 218513295 12 

f 6 = — a 2 H a 4 H a H a 8 + a H a 12 B17 

J 32 512 1024 65536 524288 8388608 v ' 

, 3 2 405 4 2625 fi 385875 8 5893965 10 160243083 , 2 

/ 8 = -o; 2 H a 4 H a 6 H a H a 10 + a 12 B18 
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Importantly, in the computation of derivatives of Laplacian coefficients with negative j values, the absolute value 
of 7 must be taken before the application of recurrence derivative formulas such as those from Eqs. (6.70)-(6.71) of 
iMurrav fc Dermottl (|1999f) . Fig. [T2l demonstrates the error incurred by truncating some of the above expansions to 
various powers of a for a ss 0.32, the value for the planets in the HD 12661 system. 

APPENDIX C 

Here we express the auxiliary variables from Eq. (|37[) as time-dependent functions of eccentricities, semimajor axis 
ratios, and terrestrial planet-giant planet mass ratio for a giant planet residing in between (case 1), exterior to (case 
2), or interior to (case 3) two terrestrial planets. For case 1: 

Q( 13 ) =MC { 2 13) a 12 (v^TI - 1) eie 3 - 2Cf 3) e 2 e 3 cos [w x - m 2 ) (CI) 

£( 13 ) = 2Cf 3) e 2 e 3 sin (roj - m 2 ) (C2) 

/ r (i2) r (23) \ 

i? (13) = 2 a 12y ^^—ej - e 2 + a 12 ^C { 2 12) ei e 2 cos (tui - tJ7 2 ) (C3) 
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For Case 3: 
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